查看原文
其他

R语言并行计算Sen+MK趋势分析

走天涯徐小洋 走天涯徐小洋地理数据科学 2022-07-17

并行计算Sen+MK趋势分析

前面的推文中讲过Sen+MK分析

上面的方法有个很大的缺点就是R语言计算过程太慢,对此我改进了前面的代码,使用terra包重写,开启并行计算,提高计算速度。

并行计算Sen+MK

  • rast读取栅格数据
  • function中长度、年份需要根据实际情况修改
  • 输出结果包含三个波段,Z值、slope和p值(详细看参考文献推文)
library(terra)
library(trend)
#输入一个文件夹内的单波段TIFF数据,在这里是历年的NDVI年最大值
flnames <- list.files(path = './ChinaYearMean/', pattern = '.tif$')
fl <- paste0("./ChinaYearMean/", flnames)
firs <- rast(fl)

#Sen+MK计算
fun_sen <- function(x){
  if(length(na.omit(x)) <34) return(c(NA, NA, NA))   #删除数据不连续含有NA的像元
  MK_estimate <- trend::sens.slope(ts(na.omit(x), start = 1982, end = 2015, frequency = 1), conf.level = 0.95) #Sen斜率估计
  slope <- MK_estimate$estimate
  MK_test <- MK_estimate$p.value
  Zs <- MK_estimate$statistic
  return(c(Zs, slope, MK_test))
}

firs_sen = app(firs, fun_sen, cores=4 )
names(firs_sen) = c("Z""slope""p-value")

#显示输出
plot(firs_sen)
writeRaster(firs_sen, filename = "./firs_sen.tif", names=firs_sen@ptr[["names"]])
计算结果

参考文献

  1. 对比R语言Raster包和Terra包栅格计算
  2. NDVI时间序列分析之Sen+MK分析全过程梳理
  3. NDVI时间序列逐像元计算Hurst指数


您可能也对以下帖子感兴趣

R 语言 SEM 笔记:结构方程模型的多组分析(Multigroup Analysis)
R语言编程秘籍:一文掌握循环结构,让你的代码更智能!
历年工企业周边年平均 PM2.5 浓度面板数据
R可视化——基于MicrobiotaProcess包进行物种差异分析!
1998~2014 年工企周边 1~5km 范围内碳排放总量面板数据

文章有问题?点此查看未经处理的缓存